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ABSTRACT 

We present a comprehensive series of dissipationless iV-body simulations to investigate the evolution 
of density distribution in equal-mass mergers between dark matter (DM) halos and multicomponent 
galaxies. The DM halo models are constructed with various asymptotic power-law indices ranging 
from steep cusps to core-like profiles and the structural properties of the galaxy models are motivated 
by the ACDM paradigm of structure formation. The adopted force resolution allows robust density 
profile estimates in the inner ~ 1% of the virial radii of the simulated systems. We demonstrate that 
the central slopes and overall shapes of the remnant density profiles are virtually identical to those of 
the initial systems suggesting that the remnants retain a remarkable memory of the density structure 
of their progenitors, despite the relaxation that accompanies merger activity. We also find that 
halo concentrations remain approximately constant through hierarchical merging involving identical 
systems and show that remnants contain significant fractions of their bound mass well beyond their 
formal virial radii. These conclusions hold for a wide variety of initial asymptotic density slopes, orbital 
energies, and encounter configurations, including sequences of consecutive merger events, simultaneous 
mergers of several systems, and mergers of halos with embedded cold baryonic components in the 
form of disks, spheroids, or both. As an immediate consequence, the net effect of gas cooling, which 
contracts and steepens the inner density profiles of DM halos, should be preserved through a period 
of dissipationless major merging. Our results imply that the characteristic universal shape of DM 
density profiles may be set early in the evolution of halos. 

Subject headings: cosmology: theory — dark matter — galaxies: halos — halos: structure — halos: 
density profiles — methods: numerical 



1. INTRODUCTION 

In the currently favored cold dark matter (CDM) 
paradigm of hierarch ical structure formation (e.g., 
iBlumenthal et ai1ll984fl . galaxy assembly initiates with 
the gravitational collapse of overdense regions into ha- 
los of dark matter (DM), the average density of which 
excee ds that of baryonic matter by roughly six to one 
(e.g.. iSnergel et ani200.1t) . Bound in the potential wells 
of DM halos, baryons cool, conden se, and form galaxies 
with a variety of prope rties fe.g.. iWhite fc Reeall97St 
IBlumenthal et afl ll984'). Understanding the nonlinear 
gravitational collapse of the DM is the necessary first step 
in gaining insight into the physical processes of galaxy 
formation. 

During the past decade the mass distribution of viri- 
alized DM halos has been the subject of extensive nu- 
merical investigations (e.g., Dubinski & Carlberg 1991; 
Navarro et al. 1996, hereafter NFW; Fukushige & Makino 
1997; Moore et al. 1999, hereafter M99; Klypin et 
al. 2001; Fukushige & Makino 2001; Power et al. 2003; 
Tasitsiomi et al. 2004; Diemand et al. 2004; Fukushige 
et al. 2004). These studies have established that the 
spherically-averaged density profiles of CDM halos can 
be described by an approximately "universal" formula 
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over the entire range of mass and length scales that these 
simulations can resolve. This universal density profile is 
comprised of an inner steep cusp with p(r) cx r~ 7 , where 
1 < 7 < 1.5, and a power-law decline p(r) cx r~ 3 in the 
outer parts. The origin of this simple form of the density 
distribution is not yet understood, but it has been con- 
jectured that merging should play a noteworthy role in 
determining DM halo profiles (e.g. 1 lW£chsl£r_e±_aU|2 002t 
IZhao et al Jl20()3 IPoeh fc Peehlesil2003 ILu et a,Ul200% 
Numerous authors have considered the effect of merg- 
ers on the evolution of the density distribution of ha- 
los a nd galaxies u s ing controlled numerical experiments 
(e.g.. IWhitellT97cl iFarouki et al.lll98l IVillumsenl 1198.1 : 
i Capelato et alJ 119951 iSyer fc White! 119981: lBarnedfL99E : 
[ Fulton fc Barnes! 120011 iBovlan-Kolchin fc Mai 12004 : 
iMoore et alJ 120041 lAceves fc Velazquez! I2005D . These 
studies employed collisionless iV-body merger simula- 
tions and reached the conclusion that the density struc- 
ture of the merger remnants is reminiscent of that of 
their progenitors. However, most of these investigations 
used idealized, non-rotating, spherical systems to repre- 
sent the halo and galaxy models and were restricted to 
a narrow range of encounter orbital energies. Further- 
more, the great majority of these earlier studies did not 
explore the effect of hierarchical merging on the evolution 
of the density distributi o n and with the notable exception 
of lAceves fc Velazquez! l)2005|) . did not address the role 
of a baryonic component and internal angular momen- 
tum in shaping the inner density profiles of merger rem- 
nants. Yet, the gravitational field in the central regions 
of galaxies is dominated by baryons and, in the standard 
paradigm of structure formation, galaxies and galaxy ha- 
los in the universe are the end result of a complex hierar- 
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chy of mergers. Elucidating the effect of internal angular 
momentum on the evolution of density profiles is funda- 
mental as previous studies have already in dicated (e.g., 
IBullock et alll2001at lAscasibar et alJl2004|) . 

Despite the inexorable increase in the dynamic range 
of current numerical simulations, several important is- 
sues remain unresolved. For example, radiative cool- 
ing and condensation of baryons in the halo cen 
ters, which is a precursor of galaxy f ormation, steep 
ens t he density distrib ution of DM iMdwjc heLal 
19801: iBarnes fc White! 119841: iBlumenthal et all 1198 ' 



Rvden fc QiinHI1987HQneriin et, 3,1.1120(141) . Tt is essential 



to ascertain whether late mergers, that take place after 
the majority of baryonic cooling has occurred, can mod- 
ify the net effect of dissipation on halo densi ty profiles, 
as con jectured bv lLoeb fc Peebles! l|2003|) and lGao et alJ 
( 2001). In addition, it is worth asking if the density dis- 
tributions of halos today are related to the density pro- 
files of their progenitors and whether halos retain any 
memory of earlier epoch s, as is expected from analyt - 
ical considerations (e.g.. iMathurf Il988[ iDehnenl 12005). 
Addressing these questions may aid in understanding 
the dependence of halo concentrat i on on halo mass 
iNavarro et alJ Il997t IBullock et all l2001bt lEke et all 
12001ft and its relation to halo mass-accretion-rates and 
formation times (Wcchslcr et al. 2002; Z hao et all 2003). 

The complexity of halo formation in a cosmological 
context, with continuous mergers, accretion, and rapidly 
changing potential wells, hinders isolation of the mech- 
anisms that drive the evolution of halo density profiles. 
We have therefore decided to examine the evolution of 
density structure in controlled mergers of DM halos and 
multicomponent galaxies using a large ensemble of dis- 
sipationless A-body simulations. We improve upon the 
work of previous authors by adopting self-consistent re- 
alizations for the DM halo models and employing realis- 
tic models for the progenitor galaxies motivated by the 
prevailing ACDM paradigm of structure formation. Our 
simulation set is carefully designed to allow an investiga- 
tion of a much larger parameter space than before and 
our primary goal is to establish conclusively the degree 
to which the density profiles of remnants retain memory 
of the structures of their progenitors over a wide range 
of inner density slopes, internal and orbital angular mo- 
menta, and in the presence of cold baryonic components. 

We conduct numerical experiments that extend and 
expand upon those of earlier studies in at least two im- 
portant aspects. For the first time, we perform a thor- 
ough study of the effect of a merger sequence on the 
density profile of DM halos in an attempt to mimic the 
hierarchical mass assembly that characterizes CDM-like 
cosmological models. Moreover, the impact of a stel- 
lar component on the density distribution of DM during 
mergers has not been addressed to date with detailed 
simulations of realistic progenitor galaxy models. As we 
illustrate below, the central density slopes of the merg- 
ing systems are notably robust and the overall shapes 
of remnant density profiles are related to a remarkable 
degree to those of their progenitors. Our results firmly 
establish that equal-mass mergers will not suffice to mod- 
ify the central density slopes and overall shape profiles 
of the DM distribution. 

The outline of this paper is as follows. In Section [21 
we present the halo and galaxy models used in this study 



and the numerical methods we employ to construct them. 
We also describe in detail our simulation techniques and 
conduct convergence tests and numerical relaxation ex- 
periments that illustrate the robustness of our results to 
numerical effects. Section contains results of the anal- 
ysis applied to binary and multiple mergers of halos and 
multicomponent galaxies. Implications and extensions 
of these findings are discussed in Section 0] Finally, in 
Sectional we summarize our main results. 

2. NUMERICAL METHODS 

2.1. Halo Models 

We consider DM halo models with den sity profiles 
that are described by the general form ijZhaol Il996t 
iKravtsov et al.lll998|) . 



p(r) 



(r/r s )T[l + (r/r s ) Q ](/ 3 -T)/ Q 



(r < r vir ) , (1) 



where 7 is the asymptotic inner slope of the profile, (3 
is the outer slope, and a parametrizes the transition be- 
tween the inner and outer profiles, with larger values of a 
corresponding to sharper transitions. Here, p s is a char- 
acteristic inner density and r s denotes the scale radius 
defined as the distance from the center where the loga- 
rithmic slope is the average of the inner and outer slopes, 
dlnp(r)/dlnr = —(7 + /3)/2. The logarithmic slope cor- 
responding to the density profile of equation JQ) is given 

by 



dlnp(r) 
d In r 



= -7-(/3- 7 ) 



(r/r s )° 
[1 + (r/r s )«] 



(2) 



It is worth emphasizing that the choice of the functional 
form of the halo profiles is not unique. Several recent 
numerical studies have suggested that the density distri- 
butions of DM halos are better represented by a function 
with a logarith mic slope continuou s ly varying as a power 
law in radius llStoehr et all 120021 iNavarro et al.ll2004t 



iMerritt et aIll20oXlGraham et all 1200.51) . However, the 
adopted profiles constitute a sufficiently realistic repre- 
sentation of the density distribution of cosmological halos 
for the purposes of this study. 

We define the virial radius, r v ; r , as the radius en- 
closing an average density equal to the virial overden- 
sity, A v ; r , times the critical density for a flat universe, 
p cr it. Thus, the virial mass and virial radius are related 
through Al™. — (4/3)7ryO cr itA v i r r^ ir . Throughout the pa- 
per we adopt the concordance ACDM cosmological model 
(Q m = 0.3, flA — 0.7, h — 0.7) and assume 2 = 0. 
The virial overdensity is then equal to A v j r ~ 103.5 
(e.g- lLacev fc Coldll993t lEke et al] fl9981. The general 
form of equation JQl includes as special cases many of 
the density profiles commonly used to fit DM halos in 
cosmological simul ations. For example, the NFW and 
iMoore et alJ l)1999|) density profiles correspond to the pa- 
rameters [a,fl, 7] = [1,3,1] and [a, /3, 7] = [1.5,3,1.5], 
respectively. 

One of the primary goals of this study is to exam- 
ine the evolution of the inner density slopes during 
equal-mass dissipationless mergers. Hence, we fix the 
outer logarithmic slope to (3 — —3, as demonstrated by 
cosmological simulation s (e.g., NFW lMoore et al1ll999l 
lAvila- Reese et aDll999|) . However, density profiles with 
outer slopes (3 > — 3 lead to cumulative mass profiles that 
diverge as r — > 00. In order to keep the total mass finite, 
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it is necessary to introduce a cutoff in the density profile. 
Truncating the profile sharply for r > r v i r , results in un- 
physical models, therefore, we implement an exponential 
cutoff which sets in at the virial radius and turns off the 
profile on a scale rjecay The truncation scale r decay is 
a free parameter and controls the sharpness of the tran- 
sition. Explicitly, we model the density profiles of halos 
beyond r V j r by 



pir) 



c(l 





T T vyc 


I exp 






''"decay 



(r > r vlr ) 



(3) 

where c = r V i r /r s is the concentration parameter. Fi- 
nally, in order to ensure a smooth transition between 
the profile interior to r v i r given by equation and the 
profile exterior to r v ; r given by equation ©, we require 
the logarithmic slope to be continuous at r = r v i r . This 
implies 

-7 - f3c a r v j r 



1 



?"dccay 



(4) 



Note that this procedure results in some additional 
bound mass beyond r v ; r , the precise amount of which 
depends upon the adopted model parameters. In the 
models we consider here, the mass exterior to r v ; r is typ- 
ically about ~ 10% of the mass contained within the 
virial radius. 

We construct iV-body hal o models for this study us - 
ing the method described in iKazantzidis et alJ (|2004b), 
which is based on an explicit computation of the exact 
phase-space distribution function (DF). Under the as- 
sumptions of spherical symmetry and an isotropic veloc- 
ity dispersion tensor, the DF depends only on the binding 
energy per unit mass £ . In this case, the DF is given by 



1 



V8ir 2 



£ d 2 p dV> 



d^y/S^P 



J_/dp\ 



(5) 

ijEddingtonl 1 1 9 1 6|) . where p(r) and ip(r) are the density 
profile and relative gravitational potential corresponding 
to the halo model, respectively. The second term on 
the right-hand side in equation (0) vanishes for any sen- 
sible behavior of ip( r ) an d p(r) at large distances. The 
d 2 p/dip 2 factor in the integrand would be difficult to deal 
with numerically, but it can be evaluated analytically us- 
ing equations Q and © f°r p{ r ) t° give an expression 
in which the only derivatives remaining are dip/dr and 
d 2 ijj/dr 2 . Both of these can be written in terms of the 
density profile p(r) and its cumulative mass distribution 

Mi?), 



d 2 P 

d^ 2 



\ GM ) 



fp 
dr 2 



dp 
dr 



Airp r 2 
M 



(G) 



Thus, equation JSJ is reduced to a simple quadrature with 
no numerical differentiation required which is integrated 
numerically to obtain the DF for the anticipated values of 
the energy 8. A Monte Carlo realization of the iV-body 
model is then readily generated by randomly sampling 
the particle positions and velocities from the DF. This 
method makes no assumptions about the functional form 
of the local velocity distribution, and therefore produces 
self-consistent equilibria that are ideal for studies of the 
detailed structure of collisionless systems. 



It is well established that cosmological halos exhibit de- 
part ures from spherical s ymmetry and velocity isotropy 
(e.g.. lCole fc Lacevlll99l| . Nevertheless, in § EH1 we an- 
alyze merger simulations between systems resulting from 
binary encounters of the initial spherical and isotropic 
models. These merger remnants are known to reproduce 
the shape and velocit y anisotropy trend s found in cos- 
mological simulations ijMoore et alJl2004|) . 

We model our halos with three density profiles specified 
by particular choices of the parameters in equation Q • 
The first model follows the NFW density profile with 
[a, (3, 7] = [1,3,1], while the second and third models 
correspond to profiles with steeper ([a,/3, 7] = [2, 3, 1.7]) 
and shallower ([a,/?, 7] = [2,3,0.2]) inner slopes. In 
the interest of brevity, throughout the remainder of the 
text, we shall refer to these density profiles as the NFW, 
"steep," and "shallow" profiles respectively. The mo- 
tivation for using the shallow density profile is three- 
fold. First, it is used to fit the inferred D M density dis- 
tribu tion in some observed systems fe.g.. Ide Blok et all 
1200 1[) . Second, it corresponds to a regime applicable 
for warm DM (WDM) and decaying DM (DDM) mod- 
els, which would introduce such a core in the profil e 
(e.g.. lHogan fc Dalcantonl2000HAvila-Reese et al.l2001(l . 
Third, it has been argued that the extrapolated density 
profile s of cosmological D M halos have cores and not 
cusps ijStoehr et al J 120021) . It is worthwhile to stress 
again that in our modeling 7 denotes the asymptotic 
slope of the density profile. Instead, the logarithmic 
slopes of the steep, NFW, and shallow profiles at 1% 
of r vir are dlnp/dlnr ~ -1.72, ~ -1.21, and ~ -0.24, 
respectively. Throughout the remainder of this paper we 
use the terms inner and central density slope to refer to 
the logarithmic slope of the spherically-averaged density 
profiles at the minimum resolved radii of the simulations. 

Our choice of density slopes ensures that the employed 
density profiles have significantly different shapes, a fact 
which will be crucial for the interpretation of the results. 
Each of our initial DM halos had a virial mass of M v - a — 
10 12 M Q , implying a virial ra dius of r v \ r — 256.7 kpc, 
and a concentration of c = 12 l|Bullock et alJl2001hfl . re- 
sulting in a scale radius of r s ~ 21.4 kpc. The adopted 
value of M v i r serves merely practical purposes and does 
not imply anything special about the particular choice of 
mass scale. Because we do not consider non-gravitational 
processes such as gaseous dissipation, the scale-free na- 
ture of gravity allows the rescaling of our models to any 
system of units and hence the extension of our conclu- 
sions to mergers between equal-mass systems of any mass 
scale. 

The left panel of Figure^shows the density, p(r), (bot- 
tom panel) and logarithmic slope profiles, G?lnp(r)/dlnr, 
(top panel) for all initial halo models as a function of ra- 
dius in units of the virial radius, r v i r . Density is given 
in units of the virial density, p V ir = A V i r p cr i t . Thick lines 
correspond to the cumulative mass profiles for all initial 
halo models, M(r), normalized to the virial mass, M v { r . 
We computed density profiles from the force resolution 
(2e, where e denotes the gravitational softening), outward 
defining the location of the most bound particle as the 
center of the system. To reduce noise, density profiles 
were averaged over 25 timesteps spanning a timescale 
equal to the crossing time at the virial radius of the sys- 
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tern, t cross (r vir ) = ^/r^JGM^. We derived logarithmic 
slope profiles by fitting the averaged density profile data 
locally about each radius. The accuracy of the density 
slope calculation was tested by comparing the results of 
the local fits against the exact analytical expressions for 
the initial models. The adopted fitting procedure yielded 
maximum relative deviations of approximately 2%. This 
accuracy suffices for the purposes of our analysis. 

2.2. Galaxy Models 

Disk galaxy models are constructed according to the 
procedure described in lHernauistl 1)199.11) and their struc- 
tural parameters are motivated by the currently fa- 
vored g alaxy formation paradigm in the standard ACDM 
model l)Mo et alJH998|) . Each disk galaxy model consists 
of a stellar disk embedded in a spherical and isotropic 
NFW DM halo. The stellar disks follow an exponen- 
tial distribution in cylindrical radius R and their vertical 
structure is modeled by isothermal sheets 

M d ( R\ 2 / z \ 

Mi? ' z) = 4^ ex H~^J sech UJ ' (7) 

where Md, i?d and Zd denote the mass, radial scale 
length, and vertical scale height of the disk, respectively. 
In our modeling, we parametrize the disk mass to be a 
fraction md of the halo virial mass, Md = mdAf v j r , and 
we specify the disk vertical scale height in units of the 
radial disk scale length. 

While most of the models used in this study do not 
contain a bulge, we also consider more general galaxy 
models that include a compact bulge component. For 
simplicity, we ignore any rotation in the bulges and model 
them as non-rotating spheroids that follow the lHernauistl 
1)1990(1 density profile 



where M b is the bulge mass and a its scale length. In 
analogy to the treatment of the disk, we assume that 
the bulge mass is a fraction nit, of the halo virial mass, 
Alb = TObM™, an d parametrize the bulge scale length in 
units of the radial disk scale length, a = /t,i?d- 

The DM halos are exponentially truncated beyond r vir 
using equation @ and have the same M v i r and c as the 
pure halo models. Once a set of cosmological parameters 
is adopted, the virial quantities M v ; r and r v ; r are uniquely 
determine d by the halo cir cular velocity at the virial ra- 
dius, Kir l)Mo et al.lll998h . Furthermore, the DM halo 
carries some net angular momentum specified by the di- 
mensionless spin parameter, A = J/Gy / \E\/M^~, where 
J and E are the total halo an gular momentum and en- 
ergy, respectively. We follow ISpringel fc White! l)1999fl 
and distribute the angular momentum of the DM halo 
by setting the halo streaming velocity to be a fixed frac- 
tion of the local total circular velocity. In our modeling, 
we assume that there is no exchange of angular momen- 
tum between the disk and the DM halo. We also assume 
no angular momentum exchange between the disk and 
the bulge during galaxy formation. Thus, the specific 
angular momentum content of the disk is conserved. 

Finally, the DM halos are adiabatically contracted in 
response to the growth of the collisionless stellar com- 
ponents under the assumptions of spherical symmetry, 



homologous contraction, circular D M particle orbits, 
and a ngular momentum conservation l)Blumenthal et all 
1986). The final DM distribution, M f (r), can then be 
derived from the initial mass profiles of DM, Mdm^), 
and baryons, Mb(r), and the final baryon distribution, 
Mb(ff), according to 

[M DM (r) + M b (r)] r = [M f (r) + M h (r f )]r f , (9) 

where is the final radius of a DM particle. The right 
panel of Figure Q displays the density and logarithmic 
slope profiles of DM after the growth of different stellar 
components. For comparison, we also show the DM dis- 
tribution of the initial NFW profile before the baryonic 
infall. Adiabatic contraction modifies the DM density 
distribution significantly. The inner profile approaches a 
power law with p(r) oc r~ 2 after including both the disk 
and bulge components. 

For all disk galaxy models considered here, we adopt 
a halo spin parameter equal to A = 0.031, which is close 
to the mean value of s pins measured in cos mological N- 
body simulations (e.g.. lBullock et al.l2001a|) . a disk mass 
fraction of md = 0.04, and a constant vertical scale height 
across the disk equal to Zd = 0.1i?d- When a bulge com- 
ponent is added, its mass and scale radius are specified 
by mb = 0.008 and /b = 0.2, respectively. The radial 
disk scale length is uniquely determined for a given 
set of parameters M v j r , c, A, md, mb, and /b- The re- 
sulting radial disk scale lengths are either i?d = 3.8 kpc 
or i?d = 3.5 kpc, depending on whether or not the 
galaxy models are bulgeless. All disk galaxy models are 
initially dynamically stable against bar formation (the 
mean Toomre Q-parameter of the disk was about 1.4) 
and their structural parameters are in accord with dy- 
namical mass mode l Al for the Milky Way presented in 
IKlvnin et all ll200l . 

Particle positions for each component are initialized 
according to the analytic density profiles. For the par- 
ticle velocities, we assume that the velocity distribution 
at any point in space is sufficiently well approximated 
by a multivariate Gaussian whose mean velocity and ve- 
locity dispersion tensor are given by the solu t ion of the 
Jeans' equations at this point (see iHernauistl l)1993|) for 
a detailed discussion of sampling velocities). A critical 
reader may note that disk galaxy models initialized under 
this scheme are not formally self-consistent. Most inter- 
esting galaxy and halo models have local self-consistent 
velocity profiles that become strongly non-Gaussian, es- 
pecially nea r their centers . Inde ed, iKazantzidis et all 
(2004c$ and ISpringel et af] l)2005|) performed numerical 
experiments of isolated, spherically-symmetric DM ha- 
los and established that the central density cusps relax 
rapidly to an inner slope significantly shallower than the 
slope t hat was initially intend e d (see IKazantzidis et alJ 
l)2004c)) and IKazantzidis et all l)2004a|) for a full discus- 
sion regarding the implications of such instabilities for 
the evolution of substructure in CDM halos and compar- 
isons of the observed satellite dynamics with the predic- 
tions of CDM models) . 

While this is strictly correct, we argue that the galaxy 
models used in this study are not seriously affected by 
the approximate scheme that we employed. The reason is 
that the growth of different stellar components forces the 
inner DM density profile to approach a power law with 
slope 7 ~ — 2. As the central density cusp becomes closer 



The Robustness of Dark Matter Density Profiles in Dissipationless Mergers 



5 




Fig. 1. — Left: Density, p(r), (bottom panel) and logarithmic slope profiles, din p(r)/dln r, (top panel) for initial halo models as a function 
of radius in units of the virial radius, r v i T . Density is normalized to the virial density, p vlr = A v j r p cr it- The solid, dashed, and dotted 
lines show results for the steep (7 = 1.7), NFW, and shallow (7 = 0.2) profiles, respectively. The initial density profiles are exponentially 
truncated beyond r v i r because they correspond to cumulative mass distributions that diverge at large radii. The density slope profiles are 
computed by fitting averaged density profiles locally at each radius. Thick lines correspond to cumulative mass profiles for initial halo 
models normalized to the virial mass, M v i T (right axis). Right: Density and logarithmic slope profiles of DM after adiabatic contraction of 
an NFW density profile in response to the growth of different stellar components. The profile before the baryonic infall corresponds to the 
solid curve. The dot-dashed, dashed, and dotted lines are the DM profile after the growth of a bulge, disk, and a bulge+disk with parameters 
equal to those of our fiducial disk galaxy model. The DM distribution responds in degrees to the size of the total stellar component included 
in the model and its central density profile approaches a power law with slope —2 after the inclusion of both disk and bulge. 



to the p(r) oc r~ 2 profile of a singular isothermal sphere, 
the local velocity distribution approaches a Gaussian. 
In this case, the accuracy with which we compute the 
velocity distribution function is significantly improved. 
In fact, here we exploit the adiabatic contraction of the 
DM to construct nearly self-consistent disk galaxy mod- 
els that do not show any discernible density evolution 
when evolved i n isolation. We note that the numerical 
experiments of iKazantzidis et aD l|20"0 4c) adopted pure 
halo models with an inner cusp of p(r ) oc r , whereas 
the simulations of lSprineel et alJ l)2005|) considered com- 
pound galaxy models with a Hernquist DM halo and no 
adiabatic contraction. 

For completeness, we also considered two-component 
elliptical galaxies comprised of a spherical non-rotating 
stellar distribution embedded in a spherical and isotropic 
NFW DM halo. The initial stellar distributions followed 
the Hernquist density profile. We constructed iV-body 
models of elliptical galaxies adopting a procedure similar 
to the one described in Section [2~71 for the realization of 
self-consistent one-component halo models. Briefly, in 
order to embed spherical stellar distributions that are in 
equilibrium inside spherical DM halos, we simply need 
to calculate the DF of component i, fi(£). This is done 
by adding the appropriate potential term to equation J5j) 



1 



V8tt 2 



£ A2 



d 2 Pi d^ 
d* 2 V£~^ 



(10) 



where pi is the density profile of component i and ^f(r) 
V'dmC?") + ^starsM is the total gravitational potential. 



Finally, the elliptical galaxy models are constructed as 
before with the particle positions and velocities initial- 
ized from the corresponding DF. Note that the density 
profile of the DM halo will be modified when the stel- 
lar distribution forms in its center and will no longer be 
described by equation We explicitly take this into 
account when calculating the halo DF. Similar to the 
disk galaxy models, we computed the change in the halo 
density structure by assuming that it responds adiabat- 
ically to the growth of the stellar component and using 
the standard form of the adiabatic contraction model of 
Blu menthal et alJ l|1986|) . 

Before adiabatic contraction, the DM halos of the el- 
liptical galaxies had the same virial mass and concen- 
tration as the pure halo models, M v i r = 10 12 M Q and 
c = 12. The mass of the stellar distribution is chosen 
to be M s tars = O.lMvir. In order to assign the remain- 
ing free parameter of the Hernquist profile, nam e ly the 
scale length a, we followed Bo vlan-Kolchin et alJ (2005) 
and used the observed relation between stellar mass and 
effective radius R c for e arly-type galaxies in the Sloan 
Digital Sky Survey data ijShen et al.ll2003|) . 

M \ 56 

lvl stars \ n 

kpc. 



R = 4.16 
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lO n M , 

For the Hernquist profile, the effective radius R c is re- 
lated to the scale length a by R c s» 1.82a, giving a fts 
2.3 kpc for our main elliptical galaxy model. In the fol- 
lowing section, we present a detailed description of all of 
the merger simulations that we undertake in this study. 
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TABLE 1 

Summary of Merger Simulations 



Initial Models 



Remnants 



Run 71 72 Orbit |dlnp/d lnr 

(1) (2) (3) (4) (5) 



M vir 
(10 12 M ) 
(6) 



(kpc) 
(7) 



d In p/d In r| 

(8) 



HRBp 


I 


I 


T-*aT"n i c 

± Oil ClUWllL 


1.16 


1.42 


288.7 


1.13 


Bpl 


x 


I 


T-*aT"n Pt/"ll 1 C 
± <J>L CIUW11L 


1.31 


1.42 


288.7 


1.30 


Bp2 


0.2 


0.2 


T-*aT"n rir^l i c 


0.29 


1.42 


288.7 


0.28 


Bp3 


1.7 


1.7 


P&r3ybolic 


1.74 


1.56 


297.7 


1.73 


hBpl 


0.2 


1 


PetrEibolic 


0.29 1.31 


1.42 


288.8 


1.23 


hBp2 


0.2 


1.7 


T-*aT"n rir\l i c 


0.29 1.74 


1.48 


292.5 


1.70 


Bri 


1 


1 


Radial 


1.31 


1.62 


301.5 


1.28 


Br2 


0.2 


0.2 


Radial 


0.29 


1.61 


301.1 


0.28 


Br3 


1.7 


1.7 


Radial 


1.74 


1.68 


305.1 


1.73 


hBrl 


0.2 


1 


Radial 


0.29, 1.31 


1.62 


301.4 


1.21 


hBr2 


0.2 


1.7 


Radial 


0.29, 1.74 


1.62 


301.6 


1.70 


Bel 


1 


1 


Circular 


1.31 


1.54 


296.5 


1.28 


Bc2 


0.2 


0.2 


Circular 


0.29 


1.57 


298.4 


0.28 


Bc3 


1.7 


1.7 


Circular 


1.74 


1.55 


297.2 


1.73 


hBcl 


0.2 


1 


Circular 


0.29, 1.31 


1.53 


296.0 


1.22 


hBc2 


0.2 


1.7 


Circular 


0.29, 1.74 


1.52 


295.2 


1.70 


Qpl 


0.2 


0.2 


Parabolic 


0.29 


3.24 


380.0 


0.27 


Qp2 


1.7 


1.7 


Parabolic 


1.74 


3.31 


382.8 


1.72 


Qrl 


0.2 


0.2 


Radial 


0.29 


3.15 


376.4 


0.27 


Qr2 


1.7 


1.7 


Radial 


1.74 


3.26 


380.8 


1.72 


Spl 


1 


1 


Parabolic 


1.31 


1.94 


320.5 


1.28 


Sp2 


1 


1 


Parabolic 


1.31 


2.61 


353.6 


1.27 


Sp3 


0.2 


0.2 


Parabolic 


0.29 


1.80 


312.1 


0.26 


Sp4 


0.2 


0.2 


Parabolic 


0.29 


2.21 


334.7 


0.25 


Sp5 


1.7 


1.7 


Parabolic 


1.74 


2.35 


341.5 


1.72 


Sp6 


1.7 


1.7 


Parabolic 


1.74 


3.09 


373.9 


1.71 


egBp 


1 


1 


Parabolic 


2.03 


1.30 


280.5 


2.01 


hbBp 


1 


1 


Parabolic 


1.57 


1.40 


287.5 


1.55 


dgBpl (coplanar) 


1 


1 


Parabolic 


1.70 


1.38 


286.1 


1.69 


dgBp2 (inclined) 


1 


1 


Parabolic 


1.70 


1.39 


286.6 


1.69 


dgBr (inclined) 


1 


1 


Radial 


1.70 


1.54 


296.7 


1.69 


hbdBp (inclined) 


1 


1 


Parabolic 


1.76 


1.38 


286.1 


1.74 



Note. — Columns (l)-(5) refer to the properties of the initial models. Columns (6)-(8) refer 
to the properties of the merger remnants. The quantities in each column are as follows. Col. 
(1): Abbreviations for the merger simulations, B = binary, Q = quadruple, S = sequential, p = 
parabolic, r = radial, c = circular, h = hybrid, eg = elliptical galaxy, dg = disk galaxy consisting 
of a halo and a disk, hb = galaxy consisting of a halo and a bulge, hbd = galaxy consisting of 
a halo, a disk, and a bulge. For mergers between galaxies containing disks the orientation of 
the merging disks is included in parentheses. All initial halo models have N = 2 x 10 5 particles 
except for run HRBp which refers to the high-resolution binary, parabolic merger between NFW 
halos with N — 2 x 10 6 particles per progenitor. Col. (2), (3): Asymptotic inner density slopes. 
Note that for galaxy models the asymptotic inner density slope of the initial uncontracted NFW 
halo model is given. Col. (4): Initial orbital configuration. Col. (5): Logarithmic density slope 
at the minimum resolved radii. For mergers of identical progenitors only one number is shown 
while for the hybrid mergers of two different initial models, two slopes are given. For the 
galaxy models we only give the density slope of the DM component. We estimate that our 
determination of slopes is accurate to within ~ 2% at radii that are well resolved. Col. (6): 
Virial mass in units of 10 12 Mq- Col. (7): Virial radius in kpc. Col. (8): Logarithmic slope 
at the innermost resolved radius. 



2.3. Description of Merger Simulations 

All numerical calculations discussed in this paper were 
carried out using PKD GRAV, a multi-stepping, paral- 
lel, tree TV-body code llStadell 120011) . PKDGRAV uses 
a spline softening length, such that the force is com- 
pletely Keplerian at twice the quoted softening length, 
and multi-stepping based on the local acceleration of par- 
ticles. We used an adaptive, kick-drift-kick leapfrog in- 



tegrator with individual particle time steps Ati chosen 
according to Ati < ^(ei/o^i) 1 ^ 2 , where e; is the gravita- 
tional softening length of the particle, ai is the value of 
the local acceleration, and 77 is a parameter that specifies 
the size of the individual timesteps and consequently the 
time accuracy of the integration. 

2.3.1. Numerical Parameters 
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For all merger simulations between pure DM halos, we 
used N = 2 x 10 5 particles and employed a gravitational 
softening length of e — 1.5 kpc. This choice enables us 
to resolve density profiles to ~ 1% of the virial radii of 
the simulated systems. The adopted force resolution is 
comparable to the mean particle separation within the 
region we want to resolve and is in accord with recent 
studies suggesting optimal scalings between the number 
of particles, N, and the minimum resolved radii (e.g., 
iPower et alll2003t iReed et al.ll2005[) . In §E3I we quan- 
tify the numerical effects of two-body relaxation on the 
central regions of our models and conduct convergence 
tests to ensure that the adopted mass and force resolu- 
tion is adequate for the purposes of this study. The re- 
sults of these tests confirm that our A^-body simulations 
are robust to artificial numerical effects. 

For binary mergers between multicomponcnt systems, 
we used N = 2 x 10 D particles to represent the DM halo 
and N = 2 x 10 4 collisionless stellar particles for each 
stellar component. Gravitational forces for the stellar 
components were computed using a softening length of 
e = 0.25 kpc and the force resolution of the DM com- 
ponent was the same as in the halo-only mergers. To 
quantify the extent of artificial numerical effects on our 
findings, we compared the results of these merger simula- 
tions against t he higher resolution calcul ations presented 
in the study of Kazantzidis et a D pOOSh . These authors 
employed galaxy models consisting of 10 6 DM particles 
and 10 5 particles in each stellar component and corre- 
spondingly smaller softening lengths. The comparison 
showed excellent agreement and in the remainder of the 
paper we will always present results for merger simula- 
tions with the standard mass and force resolution. 

For all simulations, we set the base-timestep to be 
equal to 1% of the dynamical time at the half-mass ra- 
dius of the model and allowed the individual particle 
timesteps to be at most a factor of 2 30 smaller. The 
time integration was performed with high enough accu- 
racy such that the total energy was conserved to better 
than 0.1% in all cases, which is adequate for the type of 
study that we undertake in this paper. As we have al- 
ready stated, one of our main goals is to investigate the 
evolution of the central density slopes on small scales. 
The total energy contained in the inner regions of our 
models is a few tens of a percent of that of the entire 
system, so the energy conservation accuracy must be at 
least comparable to that in order to resolve meaningfully 
the dynamics of the region of interest. 

2.3.2. Orbital Parameters 

Initial conditions for binary mergers were generated by 
building pairs of halo or galaxy models and placing them 
at a distance equal to twice their virial radii. In the 
coordinate system chosen to describe the merger simula- 
tions, the orbital plane coincides with the x — y plane, 
and the center of mass of the system coincides with the 
coordinate origin. We simulated binary mergers of sys- 
tems on parabolic orbits in agreement with co smological 
expectations (e.g.. iKhochfa r &; Burkertl l2005'). In order 
to examine the effect of encounter geometry and initial 
orbital energy on the density structure of the merger rem- 
nants we also considered circular or bits and radial orbi ts 
with zero orbital angular momenta l)Moore et aLl feOO^. 

For the parabolic mergers, we set the initial models on 



orbits with pericentric distances of 20% of the halo virial 
radii, a value that is ty pical of merging halos in cosmolog- 
ical simulations (e. g., iGhigna et aljfl998t iZentner et alJ 
120051 lBensonl l2005). The initial center of mass velocity 
of each pair was determined from the corresponding Ke- 
plerian orbit of two point masses. The trajectories of the 
merging systems deviate shortly after they begin to over- 
lap as orbital energy is dissipated and eventually follow 
closer orbits. In the case of the radial mergers, we set 
the initial relative velocities to be equal in magnitude to 
the velocities of point particles in a circular orbit about 
the common center of mass. In this way, both the radial 
and circular mergers have precisely the same total orbital 
energy. 

For binary mergers between disk galaxy models it is 
also necessary to choose a relative orientation for the disk 
components. We simulated parabolic and radial merg- 
ers using both randomly-inclined and coplanar prograde 
(disk spin vectors parallel to the orbital angular momen- 
tum vector) disk orientations. A subs et of these binary 
merger simulations were discussed in Kazantzidis et al. 
(2004aJ where additional details can be found. Finally, 
we performed a single parabolic binary merger between 
elliptical galaxies and halo+bulge systems aimed at es- 
tablishing the generality of our basic conclusions. 

Besides conducting binary mergers, one of our addi- 
tional goals was to assess the net effect of a hierarchy of 
mergers, which is a more appropriate characterization of 
the early evolution of halos in the standard paradigm of 
cosmological structure formation. Along these lines, we 
have simulated a sequence of binary, equal-mass mergers 
between identical halos for each of the shallow, NFW, 
and steep initial density profiles. For this ensemble of 
simulations we only utilized parabolic orbits. Each se- 
quence of mergers consisted of three simulations. The 
first stage was simply the binary merger between the ini- 
tial halo models described above. In the second stage, we 
merged identical copies of the merger remnants from the 
first stage. We identified a time after which the central 
density profile of the remnant did not evolve significantly 
(changes of the order of 1 — 2% were considered accept- 
able) and removed all unbound particles. We underscore 
that the out er regions of remnants ma y evolve for much 
longer (e.g., IKazantzidis et alJl2004aj) as density waves 
associated with the merger process may persist for an 
extended period of time after the cores of the halos co- 
alesce and the merger would conventionally be deemed 
complete. Finally, we chose the initial relative orien- 
tation of the principal axes of the remnants randomly, 
placed the remnants at a relative distance equal to twice 
their virial radii, and set them on parabolic orbits. For 
the third stage, we repeated the previous process using 
as progenitor systems the remnants of the previous stage. 

We determined the bound mass in ea ch of the remnants 
using t he iterative scheme described in lKazantzidis et alJ 
(2004cJ). In the rest frame of the most bound particle, 
we calculated the binding energies of all other particles 
using the tree-based gravity calculation performed by 
PKDGRAV and we removed all particles with positive 
binding energy. We repeated this calculation of binding 
energies and subsequent removal of unbound particles 
until no more unbound particles were found. In practice 
this iterative procedure converges rapidly and ensures 
that the true, bound entity will be identified. This tech- 
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Fig. 2. — Left: Spherically-averaged density, p(r), (bottom panel) and logarithmic slope profiles, d In p(r)/d lnr, (top panel) for two binary, 
parabolic merger simulations between NFW halos (runs HRBp and Bpl). Each simulation has a different force and mass resolution. In 
order to highlight any differences at small radii, we plot the product p(r)r 2 , which is constant for isothermal distributions. The density 
profiles are plotted from the force resolution (2e) outward and the number of particles, N, is indicated in the lower left-hand corner. Radius 
is shown in both physical units (bottom axis) and as a fraction of the virial radius of the initial models, r v i r (top axis). For clarity, the 
density slopes corresponding to the remnants (thick lines) have been shifted up by 0.5 dex. The downward arrow indicates the gravitational 
softening used in the low-resolution simulation. The profiles are indistinguishable from the radius corresponding to the force resolution 
of the lower resolution simulation indicating that numerical convergence has been achieved. Right: Density profile as a function of time 
spanning 18 tcross (^vir) f° r an isolated halo following the shallow density profile. In order to emphasize any differences at small radii, we plot 
the product p(r)r ( - 9+ ~'^ 2 . The maximum of this quantity indicates the scale radius of the system. Top panel: Relative density differences, 
S p = (pi n — pfl)/pi n , between the initial, pi n , and final, pg, profiles for the same timescales. Virtually no evolution in the density profile 
can be discerned over the timescales of the simulation, indicating that our numerical results are robust against two-body relaxation for 
evolution on timescales smaller than 18 tcross (r v - lr ). 



nique is essentially the same used in most group-findin g; 
routines, like the publicly-available SKID i jStadel feOOlh 
but has the advantages of using a tree structure for the 
potential computation, which requires 0(N log N) oper- 
ations instead of 0(N 2 ) for N particles, and a parallel 
implementation for very large N. In this way a much 
larger number of particles can be used in the calculation 
than would be possible with SKID at a fraction of the 
computational burden. 

In addition, we considered the simultaneous merging 
of four systems. Initial conditions for these quadruple 
mergers were generated by building halo pairs identical 
to the ones used in the binary merger experiments. We 
considered only shallow and steep initial density profiles, 
as these are the extremes of the range of profiles that we 
studied with the binary merger experiments, and set the 
center of mass position and velocity of the first pair to 
be the same as in the binary merger experiments. For 
the second pair, we changed the signs of the center of 
mass coordinates of each halo in such a way that its ori- 
entation was rotated 90° with respect to the first pair. 
We simulated quadruple mergers in which each pair was 
initially set on either a parabolic or a radial orbit. The 
trajectories of the merging systems will strongly deviate 
from a parabolic and radial orbit shortly after the begin- 
ning of the simulation, but here we are not interested in 
reproducing a particular orbital configuration. The setup 



of these experiments simply ensures that four halos that 
correspond to the steep and shallow density profiles will 
merge following orbits with very different orbital config- 
urations and orbital energies. However, for convenience, 
we shall refer to these quadruple mergers as "parabolic" 
and "radial". 

In addition to the DM-only halo mergers, we studied 
the merging of the multicomponent galaxy models de- 
scribed in § 12.21 We studied mergers of identical galaxies 
with stellar disks and bulges embedded in adiabatically 
contracted DM halos as well as mergers of galaxies with 
only stellar disks or bulges. For the halo+disk mergers, 
we studied variations in both the initial orbital config- 
uration and the initial relative orientation of the stel- 
lar disks. We considered mergers of disks randomly in- 
clined with respect to the orbital plane on both parabolic 
and radial initial orbits and mergers of coplanar pro- 
grade disks on parabolic initial orbits. For the mergers 
of halo+disk+bulge compound galaxy models the stel- 
lar disks were randomly inclined and the initial orbit 
was parabolic. Finally, for the halo+bulge and ellipti- 
cal galaxy mergers, parabolic was again the initial orbit 
of choice. Tabled contains a summary of all merger sim- 
ulations performed in this study and a list of parameters 
and variables for all initial models and remnants. 

2.4. Numerical Tests 
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In order to minimize any concern that our results 
might be compromised by artificial numerical effects, we 
performed experiments varying the mass resolution by 
a factor of 10 and scaling down the softening lengths 
according to e a jV" 1 / 3 . Results from one of these 
tests are displayed in the left panels of Figure Al- 
though the merger remnants exhibit significant depar- 
tures from spherical symmetry, it is instructive to calcu- 
late spherically-averaged remnant density profiles. Fig. [21 
shows average density p(r), and logarithmic slope pro- 
files dlnp(r)/dlnr, for two parabolic, binary mergers be- 
tween NFW halos at two different numerical resolutions 
(runs HRBp and Bpl). The agreement in the density 
and logarithmic slope profiles of the remnants simulated 
at different resolutions indicates that numerical conver- 
gence has been achieved. 

Two-body relaxation is a numerical artifact associated 
with the use of particles to sample phase-space in the 
collisionless limit. The coarse-grained sampling of phase- 
space originating from the enormous difference between 
the mass of a particle in iV-body simulations and that of 
DM particle candidates results in a mean potential that 
can be dominated by two-body interactions. As a re- 
sult, two-body relaxation sets a limit on the region of an 
./V-body system within which the numerical results can 
be trusted. All of our initial models have a central den- 
sity slope shallower than that of an isothermal (p oc r~ 2 ) 
model, therefore energy transfer due to two-body relax- 
ation would cause an expansion of the central region and 
a subsequent flattening of the inner density slope, poten- 
tially interfering with the interpretation of our results. 

To thoroughly investigate the effects of two-body relax- 
ation on our results, we performed a test simulation. We 
evolved a halo in isolation following the shallow density 
profile for a total elapsed time of t TC i ax = 18 i cross (r v i r ). 
For the particular model we studied this corresponds to 
£ ro i ax ~ 34 Gyr. The results of this experiment are shown 
in the right panels of Figure These panels present the 
density profiles, p(r), (bottom panel) and the relative den- 
sity differences, 5 P — (p m — pa)/pin, (top panel) between 
the initial p m , and final pa, profiles as a function of time. 
The density differences are small (of the order of a few 
percent), indicating that discreteness effects associated 
with the finite number of particles do not cause this par- 
ticular halo model to evolve away from its equilibrium 
configuration over these timescales. We stress that this 
shallow halo model will be the most susceptible to two- 
body relaxation owing to the relatively small number of 
particles in its central region. Therefore, we anticipate 
that all of our simulations should be unaffected by nu- 
merical relaxation for evolution on timescales that are 
smaller than t YC \ ax above. 

3. RESULTS 

3.1. Binary Halo Mergers 

We begin with the binary merger simulations between 
pure DM halos. In all parabolic and radial mergers, the 
systems merge in three to five orbits (between ~ 5 to 
~ 7 Gyr), with the specific value depending on the in- 
ner power-law indices of the systems and the initial or- 
bital configurations of the mergers. These timescales are 
larger than those report ed in earlier studies of b inary 
equal-mass mergers (e.g.. lBarnes fc Hernauistlll99ffl due 



to the larger and more realistic pericentric distances 
adopted here. The merger timescale is set by the combi- 
nation of dynamical friction, which dissipates the orbital 
energ y of the dense hal o cores, and mass loss processes 
fe.g.. lTaffoni et alJ2003lK Models with steeper inner den- 
sity distributions merge faster due to the stronger grav- 
itational drag exerted by the galaxy cores. The differ- 
ences in merging times can be of the order of ~ 1 Gyr 
or even larger between shallow and steep density pro- 
files. Finally, for circular orbits the merger process takes 
substantially longer time to complete owing to the much 
slower rate at which orbital energy is dissipated. 

The first results that we present are those of the bi- 
nary mergers involving identical halo models (runs Bpl- 
Bp3, Brl-Br3 and Bcl-Bc3). In Figure we show the 
spherically-averaged density, p(r), (bottom panels) and 
logarithmic slope profiles, d\np(r)/d\nr, (top panels) of 
the remnants for all nine identical-halo mergers along 
with the initial density profiles of each model. Note that 
each of the remnant profiles is normalized to its own virial 
radius, r v j r , as defined in § 12.11 Figure |31 demonstrates 
that both the inner power-law indices and overall shape 
of the remnant density distributions are strikingly sim- 
ilar to those of the progenitor systems. We stress that 
by the shape of the density profile, we refer to the ra- 
dial dependence of its logarithmic slope as a function of 
the scaled radius r / r V j r i a i . This basic result holds for all 
initial density profiles and orbital energies that we con- 
sidered. 

Figure 0] presents cumulative mass profiles, M(r), for 
the same set of binary merger simulations. The left and 
middle panels show mass profiles for mergers of identi- 
cal initial halo models (runs Bpl, Bp2, Brl, Br2, Bel 
and Bc2), while the right panel corresponds to the "hy- 
brid" halo mergers (runs hBpl, hBrl and hBcl) which 
we will discuss in more detail below. Clearly, the cumu- 
lative measure reveals comparably less about the details 
of the inner density profiles, but this plot does display 
an interesting fact. The merger remnants contain sig- 
nificant fractions of their bound mass, ~ 25 — 35% as 
opposed to ~ 10% for the initial models, outside of what 
we formally identify as their virial radii. Parabolic or- 
bits result in a larger fraction of particles outside r v i r 
compared to bound orbits owing to the fact that there is 
a larger amount of energy available to distribute in the 
former (see also Table ^1 . These particle fractions are 
consistent with studies of halos in cosmological simula- 
tions that find that equilibri um density profile s extend 
well beyond halo virial radii l|Prada et al.ll2005|) . 

One of the practical consequences of this re- 
sult concerns the ambiguous definition of the virial 
mass and the evolution of the virial mass in an- 
alytic models for ga l axy a nd halo formation (e.g., 
Somerville fc Primackl Il999t IZentner fc Bullock! 120031 
Benson et al.ll2003t iTavlor fc Babu1l2004t iZentner et alJ 
20051) . The results reported in Figure 0] illustrate that 
a non-negligible amount of DM in the initial halo mod- 
els is heated outside of the virial radius of the remnant. 
In fact, a few percent of the initial mass of the halos 
(between ~ 1 — 3%) becomes unbound in typical halo 
mergers. 

Analytic models typically adopt a particular definition 
of halo virial mass and then assume that the virial mass 
is strictly additive during mergers. In contrast, the mass 
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Fig. 3. — Spherically-averaged density and logarithmic slope pro- 
files for initial systems (solid lines) and remnants in binary merger 
simulations between identical DM halo models (runs Bpl-Bp3, Brl- 
Br3 and Bcl-Bc3). The profiles are plotted as a function of radius 
in units of the virial radius and each of the remnant profiles is 
normalized to its own r v i r . The asymptotic density power-law in- 
dex, 7, of the initial profiles is indicated in the upper right-hand 
corner of each bottom panel. The dot-dashed, dotted, and dashed 
lines correspond to remnants in the parabolic, radial, and circular 
initial orbital configurations, respectively. 



profiles in Figure 21 show that, for a common definition 
of the virial radius, the virial mass of the remnant is 
not simply the sum of the virial masses of the progeni- 
tors. This suggests that a more sophisticated approach 
is warranted in semi-analytic calculations that use ex- 
tended Press-Schechter merger trees to follow the mass 
assembly histories of halos. 

Following up on the findings reported in Figs.[3]and01 
are the results of the hybrid binary merger simulations 
in which we merged DM halo-only models with substan- 
tially different values of the asymptotic density power- 
law index, 7 (runs hBpl, hBp2, hBrl, hBr2 hBcl and 
hBc2). The density profiles that result from the hybrid 
merger simulations are shown in Figure |SJ This Figure 
reveals that, in equal- mass mergers between halos with 
different power-law central density indices, the remnant 
has a profile with an inner slope c lose to that of the steep- 
est of the progenitors (see also iFulton Sz Barnes! 120011 
IBovlan-Kolchin k Mai 12004(1 . Similar to the DM halo- 
only mergers of identical initial models, this result is in- 
dependent of the precise value of the initial inner power- 
law slopes and the initial orbital energy of the progenitor 
halos. In the hybrid mergers, the net result appears to 
be that the inner core of the steepest of the two initial 
systems survives essential ly intact in the remnant, an 
effect noted previously bv IBovlan-Kolchin fc Mai ( 2004) 
for mergers between NFW halos and halos with an inner 
core of constant density. 

It is interesting to stress that the isodensity contours 
in the inner region r < 0.1r V j r of a hybrid merger rem- 
nant are nearly spherical (both b/a and c/a, where a > 
b > c are the principle axis ratios, are larger than 0.8). 
Strongly prolate remnants are emblematic of parabolic 
and radial mergers sim i lar to the ones simulate d here 
Ce.g.. lMoore et al.ll2004HKazantzidis et alJl2004a|) . The 
almost spherical inner regions of the hybrid merger rem- 
nants lends support to the idea that the steepest inner 
spherical core survives, almost unaltered, in the center 
of the final system. 

In Figure El we address the relative contribution of 
the progenitor halos to the remnant density profiles in 
parabolic mergers. In the left panel we show the quan- 
tity Ap = p Tern (r)/2pi n {r), where p rem (r) and p- m (r) de- 
note the densities of remnants and initial models, re- 
spectively. In the hybrid mergers, we have defined the 
fiducial quantity 2p in (r) = pH(r) + p£(r), where p^(r) 
and (r) correspond to the initial density profiles of 
the two systems and 71 < 72. In the right panel, we 
show the relative contribution of the progenitor halos to 
the cumulative mass profiles of the remnants, AM = 
-^7i( r )/-^rem('")- Here M rem (r) is the total mass of the 
remnant and M 7l (r) denotes the mass in the remnant 
that belongs to the progenitor halo with inner power-law 
index 71 . 

As anticipated, in mergers of identical progenitor ha- 
los, the contribution to the remnant from each halo is 
equal by symmetry, but the variations in these quanti- 
ties indicate the size of variations that are to be expected 
from numerical noise. For the hybrid encounters, the 
mass from the steeper progenitor dominates the inner re- 
gions (r < 20 kpc for the [71, 72] = [0.2, 1.0] merger and 
r < 60 kpc for the [71, 72] = [0.2, 1.7] merger) of the final 
density profiles. For example, in the [71,72] = [0.2,1.0] 
merger, 80% of the mass within the inner 10 kpc of the 
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Fig. 4. — Cumulative mass profiles for remnants in binary merger simulations between halos as a function of radius in units of r v i r . All 
profiles are normalized to the virial mass of the initial halo models, M v ; r . The two leftmost panels correspond to mergers between identical 
halo models with the initial value of 7 indicated in the upper left corners (runs Bpl, Bp2, Brl, Br2, Bel and Bc2). The cumulative mass 
profiles of the initial halos are shown by the solid lines. The right panel corresponds to hybrid mergers between an NFW halo (solid line) 
and a halo following the shallow density profile (dashed line) (runs hBpl, hBrl and hBcl). The different remnants correspond to different 
initial orbital configurations. The mass distribution of merger remnants extends well beyond the virial radius (note that > 90% of mass in 
the initial models is contained within r v i r ). 




Fig. 5. — Density and logarithmic slope profiles as in Figure[3] but for binary mergers between halos with different initial central density 
power-law indices 7 (runs hBpl, hBrl, hBcl, hBp2, hBr2, and hBc2). The initial asymptotic density slopes are indicated in the upper right 
corners of the bottom panels. The steeper of the two initial profiles in each case is shown by the solid lines while the shallower of the two 
initial profiles is depicted by the dot-dashed lines. 



remnant originated from the NFW halo, while in the 
[71,72] = [0.2,1.7] merger, the steep halo contributes 
fully 90% of the remnant mass out to a distance of more 
than 20 kpc from the center of the remnant halo. Fig- 



ure also illustrates the fact that the resulting remnant 
profile is not simply a sum of the two profiles of the pro- 
genitors. The density of the remnant everywhere within 
the virial radius is smaller than twice the density of the 
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initial models. This is because, as we noted above, mass 
is dynamically heated during the relaxation accompany- 
ing the merger and subsequently moves to larger radii. 

3.2. Sequential and Multiple Halo Mergers 

Next, we consider the merger sequences (runs Spl- 
Sp6). In the sequential merger simulations, identical 
copies of the remnants from the identical halo mergers 
discussed in § 13.11 were collided with each other after 
removing the small number (~ 1 — 3% in all cases) of un- 
bound particles. We initialized the mergers on parabolic 
orbits, placed the systems at relative distances equal to 
twice their virial radii and oriented the principal axes 
of the merger remnants randomly with respect to each 
other. We repeated this process again to yield a sequence 
of three merger remnants. We shall refer to the three 
levels of mergers as "level one," "level two," and "level 
three," respectively. The level one mergers refer simply 
to the parabolic mergers discussed in S I3.1I Note that the 
pericentric distances for all merger levels are kept equal 
to 20% of the halo virial radii. This choice is motivated 
by the results of the binary merger simulations which 
demonstrated that the density structure of the remnants 
is largely insensitive to the details of the encounter or- 
bital energy. 

In all sequences of mergers, we identified a time after 
which the central density profile of the remnant did not 
evolve significantly (as before changes of the order of few 
percent were considered acceptable) . We further allowed 
the remnants to settle into equilibrium for a timescale 
equal to one crossing time at the virial radius of the sys- 
tem, i cross (rvir). Then their equilibrium state was an- 
alyzed. The density profiles and logarithmic slopes for 
two of the remnant hierarchies, starting from the shallow 
and steep initial profiles, are shown in Figure In order 
to reduce the noise in the measurement, we superposed 
the outputs from 25 snapshots between the time at which 
the merger was completed t comp , and t comp + t cmss (r vil .) . 
Not surprisingly, the remnants at all levels of the hier- 
archy exhibit faithful memories of the inner power-law 
indices of their progenitors. This indicates that, even 
in a complex hierarchy of halo mergers, the inner slopes 
of the density profiles are well preserved. Repeated dis- 
sipationless equal-mass merging cannot modify the in- 
ner power-law index or even the general shape of the 
spherically-averaged density profiles. The evolution of 
the virial mass in these simulations of hierarchical merg- 
ing is even stronger than the one portrayed in the binary 
encounters. Indeed, the fraction of particles found out- 
side of what is formally identified as frvir rises to ~ 50% 
in the third level of mergers, emphasizing the need for 
more reliable accounting in semi-analytic models. 

Following up on the evolution of the profile shapes, 
an intriguing result of the sequential merger experiments 
is reported in Figure Here, we plot the product 
r^ +7 ^ 2 /o(r), where radii and densities are normalized 
to the virial values, as a function of r/r v i r . The quantity 
ifi + 7)/2 is equal to 1.6, 2.0, and 2.35 for the shallow, 
NFW, and steep profiles respectively. There are several 
reasons to plot this quantity. First, it scales out the 
gross dependence of the halo profile on radius so that 
small changes become more apparent. Second, the radius 
at which this quantity is maximized corresponds to the 
scale radius, r s , of the density distribution. The three 



panels of Figure |S] show that the general shape of the 
remnant density profiles changes minimally throughout 
the merging sequence. In particular, the scale radius ap- 
pears to remain fixed at a nearly constant fraction of 
the virial radius of the halo at all stages. In other words, 
the concentration parameter remains approximately con- 
stant during periods of equal-mass merging. 

The strongest evolution appears to be between the ini- 
tial and the level one merger remnant profiles. How- 
ever, this shift is small and is not seen in the subse- 
quent mergers in the sequence. We attribute the differ- 
ence in the behavior between the initial model and level 
one merger remnant to the fact that the initial merger 
involves spherically-symmetric systems, while the higher 
level mergers are between triaxial halos. 

Furthermore, it is worth pointing out the presence of 
wave-like features in the density profiles at r > 0.2r v ; r . 
These features are apparent in Figure[S]as we plot a quan- 
tity that emphasizes the details of the profiles, though 
they are also noticeable in Figs. |21 El and [7| They are 
associated with "shells" of particles left over after the 
merger and do not disappear even after several billion 
years of evolution subsequent to the completion of the 
collision. Therefore, even though the mergers involve 
smooth, structureless halos, the density profiles of the 
remnants are not completely smooth. 

The last of the DM halo-only merger experiments that 
we explored involved the simultaneous collision of four 
identical halos (runs Qpl, Qp2, Qrl and Qr2). We per- 
formed these multiple merger experiments for shallow 
and steep initial profiles on both parabolic and radial 
initial orbits as described in detail in § 12.31 One might 
anticipate that in such simultaneous multiple mergers po- 
tential fluctuations are stronger and more violent, which 
may lead to a larger degree of relaxation compared to 
the binary mergers. The resulting profiles of p(r) and 
d In p(r) I d In r are shown in Figure El In this case, again, 
the merger remnants reflect the inner power-law indices 
and overall profile shapes of their progenitors and ex- 
hibit, at most, only a small shift in their scale radii rela- 
tive to their virial radii. This reinforces our basic result, 
that the density structures of the merger remnants are 
nearly the same as their progenitors. This result can- 
not be circumvented by appealing to a complex merger 
hierarchy or episodes of multiple, nearly simultaneous, 
equal-mass mergers. In all cases, the final systems ex- 
hibit a remarkable memory of the density structure of 
their progenitors. 

3.3. Galaxy Mergers 

One of the goals of this study is to assess the general 
applicability and robustness of our results and the re- 
sults of previous studies regarding the evolution of den- 
sity profiles during dissipationless equal-mass mergers. 
As a further investigation of the generality of our pri- 
mary findings, we perform numerical simulations of bi- 
nary mergers between identical disk and elliptical galax- 
ies constructed to resemble present-day systems. The 
models we adopt do not include gaseous dissipation, but 
include the various dissipationless components thought 
to dominate the mass distribution of real galaxies as well 
as adiabatic contraction of the DM halo during galaxy 
formation. These additional features allow orbital en- 
ergy to be dissipated more quickly, the merger to come 
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Fig. 6. — Left: Density fraction, Ap = p rom /2p; n , for all binary, parabolic mergers between halos as a function of radius in absolute 
units. Here p re m and p ln denote the densities of remnants and initial models, respectively. In mergers between systems with different 
asymptotic central density power-law indices, 71 and 72, we define the fiducial quantity 2p; n = pj^ + p^, where pTj and pjj? denote the 
initial density profiles of the two systems and 71 corresponds to the smaller asymptotic power-law slope. Right: Relative contribution of 
mass in the remnants of binary, parabolic mergers as a function of radius in physical units. AM is defined as AM = M 71 (r)/M Tern (r), 
where M rem (r ) is the total mass of the remnant and M 71 (r ) denotes the mass of the progenitor halo with an asymptotic power-law index 
71 in the remnant. Line types are as in the left panel. 




Fig. 7. — Density and logarithmic slope profiles as in Figure|!|] but for the sequence of binary, parabolic mergers between identical halos 
(runs Spl-Sp6). Each merger remnant was used as a progenitor in the next level of mergers after removing all unbound particles and 
randomly orienting its axis. The asymptotic density slope of the initial profiles is indicated in the upper right-hand corner of the bottom 
panels and remnant profiles are normalized to their own virial radii. The profiles from each level of the hierarchy are shown with different 
line type as indicated in the lower left-hand corner. Hierarchical merging does not serve to modify the inner slope and overall density 
profile shape down to the limit of our force resolution. 



to completion significantly sooner, and the merger rem- 



nants to have different shapes with respect to the DM 
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Fig. 8. — Density profile shapes, scale radii, and concentrations 
of the remnants in the sequence of binary mergers (runs Spl-Sp6). 
The central density slopes of the initial profiles are indicated in 
the upper left corners of each panel. In order to emphasize the 
details of the profiles and the scale radius of the transition between 
the inner and outer power laws, we plot the product p(r)r^ +7 '/ 2 
on the vertical axis rather than p(r). In each case, t he r emnant 
profiles are normalized to their own r v i r as defined in S I2.1I Under 
the assumption that each remnant can be described by the same set 
of [a, f}, 7] slopes as the initial models, which seems well justified 
in this case, the maximum of the product p(r)r(' 3 +' 1 ')/ 2 indicates 
the scale radius of the system. 

halo-only mergers ijKazantzidis et al"ll2004a|) . 



In the disk galaxy mergers (runs dgBpl, dgBp2, dgBr 
and hbdBp), two identical spirals with self-gravitating 
disks, rotating DM halos and optional compact bulges 
merge in a variety of initial orbital configurations and rel- 
ative orientations of the stellar disks. In the halo+bulge 
(run hbBp) and elliptical galaxy merger (run egBp), two 
identical systems with non-rotating DM halos and spher- 
ical stellar distributions collide on a parabolic orbit. Fig- 
ure El displays results for all galaxy mergers performed 
in this study. In the left panels of Figure we show 
the density and logarithmic slope profiles of DM in the 
halo+disk mergers with no stellar bulge component. The 
results concerning all other galaxy merger simulations are 
presented in the right panel of Figure ITU1 

The central density slopes and the overall density pro- 
file shapes of DM are maintained in mergers between 
identical systems consisting of multiple collisionlcss com- 
ponents. Hence, the presence of cold baryonic compo- 
nents does not render the chief result of this inquiry, 
that dissipationless mergers result in remnants that are 
pract ically scaled versions of the ir progenitors, ineffec- 
tual. iBovlan-Kolchin et alJ l)2005|) reached a similar con- 
clusion in dissipationless mergers of elliptical galaxies for 
several different orbital configurations and mass models. 
The overall significance of our findings lies in the fact 
that the mixing of the various collisionless components 
during the mergers occurs in a manner that preserves the 
overall density structure of the DM component. 

4. DISCUSSION 

The results presented in this paper have interesting 
implications for the formation of a universal density pro- 
file in cosmological A-body simulations. We have es- 
tablished that the characteristic power-law indices of ha- 
los are not affected by major mergers, which should be 
ubiquitous during the early rapid-mass-accretion stages 
of halo evolution. Instead, this universal profile shape 
must be something that is imprinted on halos at the 
time of the collapse of the first gravitationally-bound ob- 
jects. In hierarchical cosmological models, the amount 
of power in the initial spectrum of density perturba- 
tions is typically washed out on scales smaller than some 
critical scale. For example, in the popular case of the 
lightest supersymmetric particle this scale is set by the 
kinetic energies of the DM particles when they kineti- 
cally decouple from the th ermal fluid (e.g., iChen et alJ 
l200lHHofmann et alJl2001(l . Thus, there exists a typical 
"smallest" mass for collapsed structures, which is of or- 
der of ~ 10 -7 -10 -5 Mq for the most frequently studied 
model of supersymmetric neutralino DM. The results re- 
ported in this paper can be used to argue that these first 
smallest structures must already have density profiles 
that follow an NFW-like universal form. This conclusion 
is in agreement with findings of direct nu merical simula- 
tions. For example. iHuss et alJ l)1999f) and lShaoiro et af] 
l)2004[) find that gravitational collapse leads to NFW- 
like profiles for a w ide variety of initial condit ions, while 
iMoore et all l|1999l) and lDiemand et ail ((2005) show that 
in cosmological A-body simulations of the first virialized 
objects with masses near truncation mass of the power 
spectrum, these early halos have NFW-like density pro- 
files. 

Given that dissipationless late-time major mergers do 
not modify the shape of the profile, the problem of un- 
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Fig. 9. — Density and logarithmic slope profiles as in Figure[!|] but for the quadruple mergers between identical halos (runs Qpl, Qp2, 
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derstanding the universality of the NFW-like density 
distribution can be relegated to understanding the rea- 
son why initial, rapid gravitational collapse leads to a 
density structure that is close to NFW. However, this 
is likely to be only an approximate picture. The ac- 
tual mass assembly involves minor mergers, for which 
the density profile of the remnant will be a non-trivial 
average of the profiles of the progenitors, as smaller 
halos sink into the ce nter and mod i fy th e inner den- 
sity distribution (e.g., iSver fc White] 119981: iDekel et all 
120031 iMa fc Bovlan-Kolchinl |2004|) . Therefore, it may 
be expected that the details of merger histor ies will in- 
duce some evolution of profile shape (e.g., lOao et aD 
2005), which may explain the existence of significant 
scatte r in the inner density p rofiles of cosmological struc- 
tures ijTasitsiomi et al.ll2004l) . Nevertheless, in the early 
rapid-accretion stage of halo evolution, the fraction of 
mass contributed by minor mergers should be rela- 
tively small and the bulk of t he mass should be assem- 
bled via major mergers (e.g., ILacev fc Col^liTiM ITflfll 
iZentner fc Bullockl2003[) so the above conclusion can still 
be applicable. 

In the context of the scaling relations of DM halos in 
cosmological iV-body simulations, another intriguing re- 
sult is reported in Figure|S| This plot illustrates that halo 
concentrations remain approximately constant through 
a sequence of mergers between identical systems. The 
controlled simulations we presented do not include the 
constant redefinition of the virial radii of individual sys- 
tems that is due to the dilution of the mean background 
density with redshift. This dilution leads to a definition 
of virial radius that grows as r v i r («) oc a, where a de- 



notes the scale factor. Combining this growth of virial 
radii with our measurement of the constancy of c, leads 
to a scaling of concentration as c tx r v j r (a)/r s cx a during 
cosmological evolution. This scaling is in broad agree- 
ment with the findings of the redshift evolution of halo 
concentrations in cosmological TV-body simulations sub- 
sequent to an early phase of rapid-m ass-acc retion (e.g., 
i Bullock et alJ l2001ri IWechsler et all 12001 IZhao et alJ 
120081 iTasitsiomi et al.ll2004HDolag et al.ll2004fl. An in- 
teresting deviation from this is that IZhao et alJ (2003) 
found that during the early rapid-mass-accretion period 
of halo growth, halo concentrations appear to remain ap- 
proximately constant, suggesting qualitatively different 
behavior during the first stages of the collapse of rare 
density peaks. 

We assessed the generality of our conclusions by ini- 
tializing a sequence of parabolic encounters according to 
the conditions of cosmological mergers at redshift z — 3, 
rather than at z = 0. The initial systems followed the 
steep density profile and had the same M v i r as the stan- 
dard halo models. We scaled the concentration param- 
eter according to c cx (1 + z) -1 for fi xed mass objects 
( Bull ock et all 12001 H IWechsler et aTl I20fr3 IZhao et a.l J 
2003) and assumed that the first level of mergers ini- 
tiates at z = 3. Our modeling results in a substantially 
denser initial system than before with r v i r ~ 79.5 kpc, 
and correspondingly smaller orbital times. Indeed, the 
first level of the hierarchy concludes after only ~ 1.5 
Gyr, as opposed to ~ 5 Gyr for the z = case. The sec- 
ond and third level of mergers complete after ~ 1.8 Gyr 
and ~ 2.1 Gyr, respectively. The remnants at all levels 
of the high-z merger sequence retain excellent memories 
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Fig. 10. — Left: Density and logarithmic slope profiles of DM as in Figure [3] but for the binary merger simulations between identical 
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DM halo, a stellar disk, and a compact bulge (run hbdBp), systems comprised of a DM halo and a bulge (run hbBp), and elliptical galaxies 
(run egBp). The presence of cold baryonic components preserves the density structure of DM in mergers between identical systems. 



of the inner power-law indices and density profile shapes 
of their progenitors, in accord with the merger sequences 
of our fiducial simulations. Furthermore, the fraction 
of particles outside the virial radii of the remnants at 
each level is substantial (between ~ 30 and ~ 50%), in 
agreement with the findings reported in § 13.21 As ex- 
pected, the results of our study are insensitive to the 
overall densities and merger timescales that characterize 
mergers occ uring at d i fferen t cosmological epochs. 

Recently, ILu et all l|2005j) have argued that frequent 
major mergers associated with the rapid-mass-accretion 
phase may set the inner p a r _1 profiles of CDM halos by 
effectively isotropizing the velocity field of DM particles 
during violent relaxation. Our simulations show that the 
shape of the density profile is well preserved during equal- 
mass mergers. Violent relaxation is therefore incomplete, 
as it does not erase memory of the progenitor proper- 
ties. Further evidence to this is provided by numerical 
studies of mergers showing that the potential fluctua- 
tions necessary for violent relaxation are damped before 
the bindin g energies of particles are c ompletely random- 
ized (e.g.. lWhitefll978t iBarnes! 119921. We also observe 
that particle orbits in merger remnants, especially for 
the initial models with steep inner cusps, exhibit radial 
anisotropy at all radii. Although the initial progenitor 
models are fully isotropic (i.e., (3=1 — of/2of = 
at all radii, where a t and ay are the tangential and 
radial velocity dispersion, respectively), merger rem- 
nants exhibit anisotropy that changes from mildly ra- 



dial with ~ — 0.2 at small radii (r/r v - K < 0.1) 
to stron gly radial with ~ 0.3 — . 4 on the outskirts 
(see alsolBovlan-Kolchin fc Mall200l IMoore et alJl2004t 
IDekel et alJ2005j) . Such a radial dependence of the shape 
of is rather similar to that observed i n CDM halos 
form ed in cosmological simulations (e.g..]Cole fc Lacevl 
ll996llColm et all200ttlFaltenbacher et al.l2005jh There- 
fore, equal-mass mergers do not lead to isotropic orbits. 
Nonetheless, we have probed the effect of mergers in the 
specific regime of equal-mass encounters that take place 
over several dynamical times. The situation may be qual- 
itatively different when a large number of mergers, occur- 
ing in rapid succession, dramatically increase the depth 
of the potential ov er dynamical timescales, as assumed 
bv lLu et alJ ((2005). Such rapid evolution of potential is 
likely to be relevant for the very earliest phases of col- 
lapse. 

Our findings can be used to place constraints on the 
expe ctations for the DM dist ribut ion in galax i es and clus- 
ters. lLoeb fc Peebles! l|2003j) and lGao et"aH l)2004ft have 
put forward the hypothesis that one of the fundamen- 
tal properties of the NFW density profile is that it is 
a dynamical attractor in collisionless evolution. More 
specifically, these authors argued that if the density 
structure of a collisionless system were to deviate from 
the NFW profile, perhaps due to baryonic dissipation 
or some other process, subsequent dissipationless major 
merging would drive the density distribution toward the 
NFW form. This "attractor" hypothesis, combined with 
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observations that stars dominate gravitationally in the 
centers of galaxies, would imply shallower DM density 
profiles in the inner regions of clusters tha n predicted 
by co llisionless cosmological simulations (e.g. JSand et alJ 
2002). However, the numerical results presented in this 
study do not support the attractor hypothesis. Our sim- 
ulations show that the inner power-law indices and the 
overall shapes of density profiles remain unmodified dur- 
ing collisionless equal-mass binary mergers, even after a 
sequence of collisions, or in simultaneous multiple merg- 
ers. This conclusion persists even if the progenitor sys- 
tems contain cold stellar components with a variety of 
properties. 

An immediate consequence is that the effects of gas 
cooling and baryon condensation in halo centers, which 
steepen the inner DM density profiles, are maintained in 
subsequent dissipationl ess major merging. This is consis- 
tent with the results of lGnedin et al.l (|2004f) . who showed 
that the effects of cooling at high redshifts are retained 
in cosmological simulations to the present day and can 
be well described by the adiabatic contraction model, 
even when DM halos undergo a period of dissipationless 
evolution. 

One of the many routes to forming bright elliptical 
galaxies in hierarchical cosmological models involves dis- 
sipationless merging in their late stages of evolution. 
Though the details of this picture have yet to come 
into focus, the results reported in this study are rel- 
evant to scenarios of elliptical galaxy formation. Ob- 
served bright and faint ellipticals exhibit a dichotomy 
in their inner stellar density distributions. Bright el- 
lipticals tend to have flatter central stellar density pro- 
files compare d to faint ellipticals, which exhibit cuspy 
profiles (e.g., iFaber et al J 119971 |K ormcn ur 
findings intimate that dissipationless hierarchical merg- 
ing of faint ellipticals alone cannot alone be invoked as 
a mechanism for forming bright elliptical galaxies. Any 
such scenario would require an additional process that 
should be responsible for the flattening of the inner den- 
sity distributions. As an example, the prevalence of rel- 
atively low-density cores in bright elliptical galaxies may 
be a consequence of the existence of supermassive black 
holes and their interactions with the stellar backg rounds 
UMakino fc FhisuzakillT^lMerntt fc CnidEnMji. 

An analytical framework exists that aids in the un- 
derstanding of our numerical results. Undoubtedly, our 
dissipationless simulations must conform to Liouvillc's 
theorem demanding that the phase-space density of a 
dynamical system that undergoes collisionless evolution 
to a new state cann ot increase. Nevertheless, as h as also 
been emphasized by Bovlan-Kol chin fc Mai l|2004|) . Liou- 
ville's theorem applies to a six-dimensional quantity and 
there is no a priori reason to expect that both configu- 
ration and momentum space distributions should be in- 
dependently conserved. In general the distributions in 
these subspaces are not conserved. In fact, our simula- 
tions show that both configuration and momentum space 
distributions are altered during mergers. 

More complicated treatment of phase-space evolution 
is needed to draw conclusions about evo lution of den- 
sity _profiles (IMathurl Il98l IDehnenl l2005j) . In particu- 
lar, |Dehne^ (|2005T) recently showed that for phase-space 
density /, the excess mass function D(f), defined as 
the difference between the mass with phase-space den- 



sity > / and the product of / and the volume of phase- 
space with density > /, can only decrease upon mixing. 
This mixing theorem and the simple properties of D(f) 
as / — > oo for self-gravitating density cusps can then be 
used to prove that a merger remnant cannot have a den- 
sity cus p steep er than that of the steepest progenitor. 
iDehnenl ((2005) also argued that reducing the slope of 
the density cusp would require arbitrarily large dilution 
of the phase-space density so that it seems implausible 
that the remnant could have an inner cusp shallower than 
the steepest progenitor cusp. 

Our findings are in agreement with the conclusions of 
IDehnenl ([2005). Indeed, Figure illustrates that in a 
merger of two identical power-law density profiles, the 
remnant exhibits a power-law inner density slope equal 
to that of the progenitors. In addition, the results em- 
bodied by Figure[S]hint that in mergers between halos of 
very different inner density power laws, the remnant has 
a central density slope slightly shallower than that of the 
steepest of the progenitors on scales that are o f practi- 
cal co ncern for galaxies. However, the results of lDehnenl 
( 2005) apply asymptotically as r — > and our simula- 
tions simply may not resolve sufficiently small radii for 
the asymptotic result to hold. In this regard, we see no 
evidence that the second conclusion of lDehnenl l)2005|) is 
violated but higher-resolution numerical simulations of 
merging density cusps are required to test this conjec- 
ture more extensively. 

Our findings also shed light on the r esilience of density 
cusps t o gravitational interactions. iKazantzidis et all 
(2004c) elucidated the dynamical effect of tides on the 
central density profiles of cuspy satellite halos. These 
authors employed collisionless, high- resolution A-body 
cosmological simulations with simulations of the tidal 
stripping of individual satellites orbiting within a static 
host potential and showed that the density distribution 
retains its initial inner power-law index as the satellite 
experiences mass loss. Their analysis revealed that the 
central density cusp is notably stable to tidal shocks and 
gravitational stripping. The results reported in this pa - 
per combined with those of Kazantzidis et al. (2004c), 
suggest that density cusps, once formed, are extremely 
difficult to disrupt or even to modify in any significant 
way. 

The preservation of the inner power laws of halo 
density profiles has intriguing implications for tests 
of the nature of DM through observations of galactic 
and sub-galactic structure. I n alternatives to CDM 
models such as WDM Ce.g.. iPagels fc Primackl 119821: 
IColombi et al.lll996t IHogan fc Dalcantonll2000l) i or DDM 
fe.g.. lKaplinghatll2005t lCembranos et alJl20*05|) the DM 
particles typically have non-negligible momenta which 
leads to a suppression of small-scale structure for two 
reasons. The first reason is that the power spec- 
trum of density fluctuations on small scales is damped 
by the free-streaming of the part icles (e.g.. iMal Il996t 
IBode et aJJl2^IKaplinghal2005|l .The second reason is 
that thes e particles are restricted to finite phase-space 
densi t ies llTremaine fe Gunnlll979t iHogan fc Dalcantonl 
120001: iKaplinghatl 12005^ The phase-space constraint 
is thought to induce a cored inner profile of nearly 
constant density in small and early-forming DM halos, 
because the very high phase-space densities of cusp- 
like profiles are unachievable (e.g., lAvila-Rccse et all 
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|2001[) . The findings of the present study imply that 
if the small, early-forming building blocks of larger ha- 
los were formed with cored density profiles, these cores 
should be well-preserved during the hierarchical assem- 
bly of larger halos, in the very least during periods of 
rapid major mergers. The signatures of the alternative 
DM could, therefore, be possibly observed in accurate 
measurements of the rotation cur ves of DM-dom i nated 
galaxies at small scales (e.g.. IdeBlok et al.l 120011; 
Salucci et all l2003t iWeldrake et all I200& iSimon et all 
2005|) or the flux anomalies in strong-lens systems 
Ce.g.. iDalal fc Kochanekl 120021: iZentner fc Bullock! 120031: 
iRozo et alTl2005| ) . This could potentially provide funda- 
mental constraints on the properties of the DM particles. 

We conclude by emphasizing that in this article, we 
have only addressed the evolution of density profiles in 
mergers between equal-mass systems. Recent works have 
examined more general situations involving systems with 
a spectrum of masses and reached similar conclusion s 
l|Bovlan-Kolchin fc Mall2004t lAceves fc Velazauezl^OO.^ . 
This fact allows us to argue that our main findings can 
be generalized to all unequal-mass major mergers with 
some degree of confidence. 

In our modeling, we have also implicitly assumed that 
cuspy halos have a higher phase-space density in the in- 
ner regions than cored profiles. The fact that in hybrid 
mergers the inner density slope of the remnant is closer 
to that of the steepest of the progenitors suggests that 
the quantity that drives the evolution of the profile is 
the phase-space density. However, it is not difficult to 
imagine merging scenarios where this situation can be re- 
versed as, for example, in cases of unequal-mass mergers 
or encounters between systems with considerably differ- 
ent concentration parameters. We defer a detailed nu- 
merical study of these considerations to future work. 

5. SUMMARY 

Using a large suite of controlled, dissipationless N- 
body simulations we have explored the evolution of DM 
density profiles in equal-mass mergers of DM halos and 
multicomponent galaxies. A common feature of the 
majority of previous related investigations is that they 
did not consider the effect of hierarchical merging that 
characterizes CDM models, nor did they examine the 
role of baryonic components on the density structure 
of merger remnants. In contrast, our numerical exper- 
iments are performed in degrees of increasing complex- 
ity. We present an ensemble of simulations comprised not 
only of binary encounters, but also sequences of mergers 
and collisions between multiple progenitors in order to 
study the range of behaviors that are realized in the con- 
text of cosmological structure formation. We have also 
probed a wider range of parameter space with respect 
to earlier studies by varying the initial orbital configu- 
rations and orbital energies in an effort to establish the 
generality of our main results. Finally, we examined col- 
lisions between systems with multiple components con- 
sisting of rotating DM halos, stellar disks, and stellar 
spheroids. This enabled us to elucidate the impact of 
internal angular momentum and the presence of cold 
baryons on our findings. We find that the merger rem- 
nants in our simulations conform to a set of several simple 
and general results that can be summarized as follows. 



1. Binary mergers between identical DM halos with 
density structures described by various asymptotic 
power-law indices par 1 ranging from steep cusps 
to core-like profiles produce remnants with inner 
density slopes equal to those of the progenitors. 
Furthermore, the overall shape of the remnant den- 
sity distribution constitutes a remarkable reflection 
of that of the initial systems. If the progenitor halos 
are constructed with appreciably different asymp- 
totic power-law indices, the inner density slope of 
the remnant is closer to that of the steepest of the 
initial systems. These conclusions hold for a wide 
range of orbital energies. 

2. The aforementioned findings remain valid for a 
variety of encounter configurations including se- 
quences of several consecutive merger events, de- 
signed to mimic hierarchical merging, and simul- 
taneous collisions of several systems. Overall, the 
inner slopes and shapes of density profiles are re- 
markably robust during dissipationless evolution, 
regardless of the number of mergers and initial con- 
ditions associated with the encounters. Dissipa- 
tionless equal-mass mergers do not provide a mech- 
anism for driving the density profiles of collisionless 
matter toward a universal form, contradicting the 
hypothesis that t he NFW profile behaves as a dy- 
nami cal attractor l)Loeb fc Peeblesl2008HGao et alJ 

Hoi- 

3. In binary mergers between identical systems, par- 
ticles from both systems contribute equally to the 
final density profile at all radii, as must be the 
case due to symmetry. In mergers between systems 
of markedly different inner density power laws, 
the central regions of the remnants are dominated 
by particles originating from the system with the 
steeper inner density slope. The steepest inner core 
survives, almost unaffected, in the remnant. 

4. Remnants in mergers between equal-mass systems 
contain significant fractions of their bound mass, 
~ 25 — 50%, well beyond their formal virial radii 
and extending out to sa 2 — 47- v i r . This suggests 
that the virial mass is not simply additive in merg- 
ers, as is commonly assumed in many semi-analytic 
models of halo and galaxy evolution. Semi-analytic 
models may significantly overestimate the mass 
and density within the virial radii of merger rem- 
nants. 

5. Mergers between identical DM halos produce self- 
similar remnants in the sense that their scale radii, 
r s , defined as the radii at which the density pro- 
files transition from an inner to an outer power law, 
are very similar to those of the progenitors when 
measured in units of the halo virial radii, r v ; r . In 
other words, halo concentrations remain essentially 
unmodified during equal-mass encounters and the 
remnants correspond to scaled versions of their pro- 
genitors. 

6. Mergers between identical systems containing cold 
stellar components in the form of disks, spheroids, 
or a combination of the two, in addition to the 
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extended DM halos, exhibit a degree of self- 
similarity akin to those of the DM halo-only merg- 
ers. Namely, the overall shape of the density pro- 
files and the inner density slopes of DM are main- 
tained in the remnants. Particle mixing preserves 
the density structure of the DM component. The 
effects of gaseous dissipation on the DM profiles of 
the progenitor systems will be retained in the den- 
sity distribution of their descendants under dissi- 
pationless major merging. 
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